Handling the peculiarities of marine movement data

1 Credits

This tutorial uses materials that was created by the CANSSI Collaborative Research Team led by Vianey Leos Barajas and Marie Auger-Méthé. I thank Fanny Dupont, Ron Togunov, Natasha Klappstein, Arturo Esquivel, Marco Gallegos Herrada, and Sofia Ruiz Suarez for their contribution to this original material.

2 Tutorial goals

The goal of this tutorial is to explore ways to pre-prepare marine movement data. The primary learning objectives are to:

1. Regularize movement tracks that are irregular using Fastloc GPS data as an example.
2. Create predicted tracks from error-prone data using Argos data as an example.

3 General setup

First, let’s load the packages that we will need to complete the analyses. Of course you need to have them installed first.

library(momentuHMM) # Package for fitting HMMs, we use it for crawlWrap
library(tidyverse)  # data management
library(ggspatial)  # plot the data
library(sf)         # spatial data processing
library(kableExtra) # produce visually appealing tables
library(geosphere)  # to calculate step lengths from lat/lon locations - just need it installed
library(conicfit)   # needs to be installed for crawlWrap
library(here)       # To help with sourcing
library(adehabitatLT) # setNA
library(terra)
library(tidyterra)
library(aniMotum)
library(mitools)  # used by mometuHMM, just need to be installed

# install.packages("aniMotum",
#                  repos = c("https://cloud.r-project.org",
#                  "https://ianjonsen.r-universe.dev"),
#                  dependencies = TRUE)

We are assuming that you are working from Madeira-Workshop repository main folder, and using here() to help find the good directory for each files. We will need the functions in the following file utility_functions.R.

source(here("D1-data-prep-ssm",
            "utility_functions.R"))

4 Regularising Fastloc GPS data

4.1 Narwhal movement data

We will analyze a dataset containing movement tracks of three narwhals tagged with Fastloc-GPS tags. The dataset was provided by Dr. Marianne Marcoux (Fisheries and Oceans, Canada). Dr. Marcoux provided the data only for this tutorial, please do not use for other purposes without their consent. Contact: Marianne.Marcoux@dfo-mpo.gc.ca.

For simplicity, we only examine the fastloc-GPS data from one week in August 2017.

Photo by Paul Nicklen

First, let’s import the raw Fastloc GPS narwhal data and convert the time column to an appropriate date format.

tracks_gps_O <- read.csv(here("D1-data-prep-ssm", "data", "tracks_fastloc_gps.csv")) %>%
  mutate(time = ymd_hms(time),
         ID = factor(ID))

We can have a first look at the data, to get more familiar with it.

head(tracks_gps_O)

##        ID                time        x       y loc_class
## 1 T172062 2017-08-08 00:00:00 -80.6414 72.2694       GPS
## 2 T172062 2017-08-08 00:20:00 -80.6525 72.2728       GPS
## 3 T172062 2017-08-08 00:42:00 -80.6665 72.2719       GPS
## 4 T172062 2017-08-08 00:52:00 -80.6692 72.2674       GPS
## 5 T172062 2017-08-08 01:09:00 -80.6682 72.2582       GPS
## 6 T172062 2017-08-08 01:19:00 -80.6613 72.2574       GPS

The columns/variables are: - ID: Individual identifier - time: time of location - x: longitude - y: latiture - loc_class: all GPS

The data we obtain is often quite messy with some rows/records missing information and other records duplicated. We can filter records to keep only complete location data using !is.na(x) & !is.na(y). To remove duplicate records (same time, location, and location class), we will use the lag function from dplyr, which will use the value from the previous row so that we can compare to the current row.

tracks_gps_O <- tracks_gps_O %>% 
  # remove missing locations
  filter(!is.na(x) & !is.na(y),
         # remove identical records
         !(time == lag(time) & 
             x == lag(x) & 
             y == lag(y) & 
             loc_class == lag(loc_class)))

Let’s plot the data over the bathymetry (i.e., depth of the ocean) and land layers, to get a sense of where these narwhals are swimming.

land <- st_read(here("D1-data-prep-ssm", "data", "NorthBaffin.shp"))

The raw bathymetry data used to create this raster was taken from the the GEBCO global ocean and land terrain model from https://pressbooks.bccampus.ca/ewemodel/chapter/getting-bathymetry/. Note that land values of the bathymetry layer are 0s.

bathy <- rast(here("D1-data-prep-ssm", "data", "bathy_4_narwhals.tiff"))

Let’s plot the narwhal movement data over the bathymetry and land layers.

ggplot() +
  geom_spatraster(data = bathy) +
  geom_sf(data = land, fill = "beige") +
  geom_spatial_path(data = tracks_gps_O, 
                    aes(x = x, y = y, colour = ID), crs = 4326) +
  coord_sf(expand = FALSE)

As we will see during the week, many analyses require that the locations are at regular time intervals, however Fastloc GPS data is often taken at irregular time intervals, since it depends on when the animal surface. When the data is irregular, there are two key decisions we must make, (1) the temporal resolution to use, and (2) how to address large data gaps.

4.1.1 Selecting a time interval (resolution)

The desired resolution depends on the biological question you are asking, as different behaviours and biological processes occur at different spatial and temporal scales (e.g., seasonal migration, daily movement between foraging and resting grounds, and fine scale foraging decisions), and the resolution of the raw data you have. Generally, higher resolution data is preferred as it has more information. However, it is possible to have too-high of a resolution wherein information from fine-scale variability drowns out the signal from coarse-scale patterns of interest (e.g., seasonal migration). A very coarse rule of thumb, is that you want 3-50 data points per behaviour. For behaviours spanning several hours, that roughly corresponds to a desired resolution between 2 min and 60 min.

First, let’s calculate the time difference between successive records using difftime and lead (compares current row to following row) and place these values in a new column called dt. Note that the time difference is in minutes (units = "mins"). For the last record of each individual (i.e., when ID != lead(ID)), we will set the value to NA.

# Calculate time difference between locations
tracks_gps_O <- tracks_gps_O %>%
  mutate(dt = ifelse(ID == lead(ID), # If next data row is same individual
                     # calculate time difference
                     difftime(lead(time), time, units = "mins"), 
                     NA))

Let’s see what resolutions may be possible in the data by looking at the most frequent time gaps.

# Visualise time differences (all and zoomed)
par(mfrow = c(1, 2))
hist(tracks_gps_O$dt, 1000, main = NA, xlab = "Time difference (min)")
hist(tracks_gps_O$dt, 1000, main = NA, xlab = "Time difference (min)",
     xlim = c(0,100))

# identify the most frequent dt
tracks_gps_O %>% 
  {table(.$dt)} %>% 
  sort(decreasing = TRUE) %>% 
  head()
## 
##  10  11  12  22   9  13 
## 400 145  90  64  63  63

We see that the most frequent time gap is \(10\) min, followed by \(11\), \(12\), \(22\), \(9\) and \(13\) min. We also see the majority of the gaps are \(< 60\) min, however some are in excess of \(600\) min. Finer resolutions will contain more data gaps. For some analyses, frequent and large data gaps can be difficult to handle, especially as the number of missing data points approaches or exceeds the existing data. Let’s examine the potential data structure at different resolutions for the different animals.

We can now use the function p_na (in the script utility_functions.R) to look at the proportion of NAs we would get with 10, 20, 30, and 60 min resolution.

# summarise track dt
tracks_gps_O %>% 
  group_by(ID) %>% 
  summarise(p_NA_10m = p_na(time, 10),     # 10 min 
            p_NA_20m = p_na(time, 20),     # 20 min 
            p_NA_30m = p_na(time, 30),     # 30 min 
            p_NA_60m = p_na(time, 60)) %>% # 60 min
  # return formatted table
  kable(digits = 3, col.names = c("Narwhal ID", paste(c(10,20,30,60), "m"))) %>%
  kable_styling() %>% 
  add_header_above(c("", "Resolution" = 4))

	Resolution
Narwhal ID	10 m	20 m	30 m	60 m
T172062	0.486	0.236	0.191	0.152
T172064	0.502	0.203	0.120	0.074
T172066	0.488	0.230	0.185	0.160

Here we see that the \(10\) min interval, around \(50\%\) of the locations would be missing.

This is a lot! In many cases, we may want to be more conservative and use a \(30\) min resolution or \(60\) min resolution.

4.1.2 Handling data gaps

There are several ways to deal with data gaps:

1. Split tracks
2. Interpolate locations
3. Fill the gaps with NAs
4. Multiple imputation

The method to use will depend on the main analysis you are interested in.

4.1.3 Splitting tracks

One way to account for missing data is to split the track where there are large gaps (i.e., assign each track segment a new individual ID). This strategy is particularly appropriate when you have long enough gaps for which interpolation method are unlikely to perform well. We can split the tracks when the gaps larger than a predetermined threshold.

Here, we will use a function (found in utility_functions.R) to split the tracks. We define the maximum allowable gap (at which point it will start a new segment), as well as the shortest allowable track segment.

These are somewhat arbitrary decisions, and depend on your subsequent choices for regularisation. In this tutorial, we will be interpolating missing locations (within each segment) and so we only want to allow gaps that can reasonably be predicted.

We are using a 30 min resolution, so we allow a 90 minute gap (i.e., we assume we can predict 2 missing locations), and we want each segment to be at least 120 min (i.e., have at least 5 locations) long so that we have enough information about state transitions.

# Use function from utility_function.R to split data at gaps
track_split <- split_at_gap(data = tracks_gps_O, 
                           max_gap = 90, 
                           shortest_track = 120)

The new data has an updated ID column, where with the original ID and track segment number. The original ID is now in ID_old

head(track_split)

##          ID                time        x       y loc_class dt  ID_old
## 1 T172062-1 2017-08-08 00:20:00 -80.6525 72.2728       GPS 22 T172062
## 2 T172062-1 2017-08-08 00:42:00 -80.6665 72.2719       GPS 10 T172062
## 3 T172062-1 2017-08-08 00:52:00 -80.6692 72.2674       GPS 17 T172062
## 4 T172062-1 2017-08-08 01:09:00 -80.6682 72.2582       GPS 10 T172062
## 5 T172062-1 2017-08-08 01:19:00 -80.6613 72.2574       GPS 10 T172062
## 6 T172062-1 2017-08-08 01:29:00 -80.6592 72.2563       GPS 15 T172062

This data is still irregular, but now we have smaller segments split when there are large data gaps. Let’s visualize the different segment.

ggplot() +
  geom_spatraster(data = bathy) +
  geom_sf(data = land, fill = "beige") +
  geom_spatial_path(data = track_split, 
                    aes(x = x, y = y, colour = ID), crs = 4326) +
  coord_sf(expand = FALSE)

Splitting the tracks is often be a first step, before interpolating or other adjustments.

4.1.4 Interpolation (correlated random walk)

Once the track is split, there is often still irregularity within each segments, and we want to interpolate or predict new locations to form a regular grid of observations.

The simplest approach is to use linear interpolation between missing times, but a better option is to predict locations from a continuous-time correlated random walk (CTCRW). momentuHMM contains wrapper functions to interpolate missing locations by fitting a CTCRW to the data based on the crawl package by Devin Johnson and Josh London. There are many options to fit the CTCRW model, and a detailed tutorial for analysis with crawl is available here: https://jmlondon.github.io/crawl-workshop/crawl-practical.html. Let’s try to fit the most basic model using the wrapper function crawlWrap. In the most basic form, we only need to provide tracking data with the columns ID, time, x, and y and specify the desired temporal resolution.

First, let us transform the data into an sf object. crawlWrap can also take a data.frame as an argument but that would imply renaming some of our columns. It is easier to just transform the data into an sf object.

Here it’s important to project the data in a good coordinate system, since the model takes into account speed/step length.

track_split_sf <- track_split %>%
  st_as_sf(coords = c("x", "y")) %>% # converts to an sf object
  st_set_crs(4326) %>% # define CRS
  st_transform(2962) # reproject data to a UTM

Now we can fit the CTCRW to each track segment and create a data frame of predicted locations. We decided above that a good, maybe conservative time interval was 30 min. We will use that interval here. Since interpolating when they are large data gaps can bias our analysis, we use the split tracks.

# crawl can fail to fit periodically, so I recommend always setting a seed 
set.seed(12)

# fit crawl
crwOut <- crawlWrap(obsData = track_split_sf, timeStep = "30 mins", theta = c(7, 0))
## Fitting 20 track(s) using crawl::crwMLE...
## Individual T172062-1...
## Individual T172062-2...
## Individual T172062-3...
## Individual T172062-4...
## Individual T172062-5...
## Individual T172062-6...
## Individual T172062-7...
## Individual T172064-1...
## Individual T172064-2...
## Individual T172064-3...
## Individual T172064-4...
## Individual T172064-5...
## Individual T172064-6...
## Individual T172066-1...
## Individual T172066-2...
## Individual T172066-3...
## Individual T172066-4...
## Individual T172066-5...
## Individual T172066-6...
## Individual T172066-7...
## DONE
## 
## Predicting locations (and uncertainty) at 30 mins time steps for 20 track(s) using crawl::crwPredict...
## DONE

Let’s look at a few track segments with interpolated values.

plot(crwOut, animals = "T172062-3", ask = FALSE)

plot(crwOut, animals = "T172064-5", ask = FALSE)

plot(crwOut, animals = "T172066-4", ask = FALSE)

Notice how the predicted tracks do not make perfectly straight lines through missing sections, which is an improvement on what a simple linear interpolation method would provide.

We can now extract the predicted regular tracks.

# Get predicted tracks from crawl output
track_int <- crwOut$crwPredict[which(crwOut$crwPredict$locType == "p"),
                              c( "ID", "mu.x", "mu.y", "time")]
colnames(track_int) <- c( "ID","x", "y", "time")
head(track_int)
##           ID         x       y                time
## 1  T172062-1 -285113.8 8176359 2017-08-08 00:20:00
## 4  T172062-1 -285813.8 8176130 2017-08-08 00:50:00
## 8  T172062-1 -286046.5 8174867 2017-08-08 01:20:00
## 11 T172062-1 -286045.6 8174408 2017-08-08 01:50:00
## 14 T172062-1 -285394.1 8174699 2017-08-08 02:20:00
## 17 T172062-1 -284706.5 8173150 2017-08-08 02:50:00

Note here that the time is at every 30 min.

4.1.5 Place NAs in (small) gaps

Instead of interpolating, one can leave data gaps as NAs. Some analysis, for example hidden Markov models without spatial covariates extracted from the locations, could easily handle NAs in the locations. The maximum size of a gap to allow depends on the frequency of the missing data, frequency of locations, study species, and behaviours of interest. However, for many analyses that assume regular time interval, it’s important to explicitly include the NAs, so that the function knows that there are missing data they have to account for. Placing NAs for missing locations can be used on its own or in conjunction with splitting tracks. The package adehabitatLR has a function setNA dedicated to it.

We will apply this to the split tracks that we used previously used in the tutorial. Here, instead of using crawl to interpolate missing locations, we are simply creating a dataframe with the missing locations set to NA (i.e., creating a regular grid of observations with some NAs). The idea is that we still want to separate tracks into segment when there are large data gaps, but place NAs when it’s just a few locations missing.

The first step is to create an adhabitat trajectory. Here, we simply input coordinates, time, and the ID (here track segment id).

track_ltraj <- as.ltraj(xy = track_split[, c("x", "y")], 
                                   date = track_split$time, 
                                   id = track_split$ID)

Now we set the NAs. As before, we use a time interval of 30 min. We specify that we have a tolerance for the imprecision in the data, here at most 15 min. We use the first location at the time reference.

# Create adehabitat object, containing the trajectory padded with NAs
track_NA <- setNA(ltraj = track_ltraj, 
                  date.ref = track_split$time[1], 
                  dt = 30, tol = 15, units = "min")

# Transform back to dataframe
track_NA <- ld(track_NA)[, c("id", "x", "y", "date")]
colnames(track_NA) <- c("ID", "x", "y", "time")
head(track_NA)

##          ID        x       y                time
## 1 T172062-1 -80.6525 72.2728 2017-08-08 00:20:00
## 2 T172062-1 -80.6692 72.2674 2017-08-08 00:52:00
## 3 T172062-1 -80.6592 72.2563 2017-08-08 01:29:00
## 4 T172062-1 -80.6558 72.2550 2017-08-08 02:04:00
## 5 T172062-1 -80.6316 72.2582 2017-08-08 02:26:00
## 6 T172062-1 -80.5892 72.2403 2017-08-08 03:02:00

We can see now that there are NAs for some locations, and that the time is not exactly at every 30 min, but close to. So this may not be adequate for some analysis that require strict time regularity. One could round to the nearest 30min, but this could be problematic.

4.1.6 Multiple Imputation

For many analyses where you need to extract covariates from the location data, it may be wise to use multiple imputation, where you create multiple tracks with a CTCRW model and extract the covariates for each of the tracks. Then you can propagate the uncertainty in other analyses.

Multiple imputation works by first fitting a CTCRW model to the original data, second, drawing (i.e., simulating) a number of realisations of the position process based on the CTCRW model, third (not done here), fitting your model of interest to each of the simulated realisations, and finally, pooling the estimated parameters. momentuHMM has several functions to implement multiple imputation. The function MIfitHMM can be used both to simulate realisations of a fitted CTCRW and fit hidden Markov model (HMMs) to each one. The number of simulations is specified with nSims. We can simulate realisations without fitting HMMs by setting fit = FALSE.

Here, let’s use first fit a CTCRW model on segmented tracks created above. We will simulate 4 tracks using MIfitHMM.

set.seed(12)

# Fit the correlated random walk, MIfitHMM takes a crwData object
crw_gps_30 <- crawlWrap(obsData = track_split_sf, timeStep = "30 mins")

## Fitting 20 track(s) using crawl::crwMLE...

## Individual T172062-1...

## Individual T172062-2...

## Individual T172062-3...

## Individual T172062-4...

## Individual T172062-5...

## Individual T172062-6...

## Individual T172062-7...

## Individual T172064-1...

## Individual T172064-2...

## Individual T172064-3...

## Individual T172064-4...

## Individual T172064-5...

## Individual T172064-6...

## Individual T172066-1...

## Individual T172066-2...

## Individual T172066-3...

## Individual T172066-4...

## Individual T172066-5...

## Individual T172066-6...

## Individual T172066-7...

## DONE

## Predicting locations (and uncertainty) at 30 mins time steps for 16 track(s) using crawl::crwPredict...

## DONE

# simulate 4 realisations of the 30 min GPS CTCRW model
MI_sim_gps <- MIfitHMM(crw_gps_30, nSims = 4, fit = FALSE)

## Drawing 4 realizations from the position process using crawl...

## DONE

## DONE

This will return warning messages. They are not shown here, but one should look at the segments that resulted in warning messages. This can potentially be fixed by forcing longer track segments above.

Let’s look at one segment: T172066-6.

ID_seg <- "T172066-6"
# plot locations for first narwhal
# filter first ID from original data
track_eg <- track_split_sf %>% 
  mutate(x = st_coordinates(track_split_sf)[,"X"], 
         y = st_coordinates(track_split_sf)[,"Y"]) %>% 
  filter(ID == ID_seg)

# filter first ID for each simulation
sim_tracks <- lapply(MI_sim_gps$miData, function(x){
  filter(x, ID == ID_seg)})

# plot original track for first narwhal
plot(track_eg$x, track_eg$y, col = "red", xlab = "x", ylab = "y", asp = 1, type = "l")

# plot each simulated track
mute <- mapply(function(data, col) {
                 points(y ~ x, data = data, col = col, type = "l")
               }, data = sim_tracks, 
               col = list("cadetblue1", "cadetblue2", "cadetblue3", "cadetblue4"), 
               SIMPLIFY = FALSE)

Notice how in some areas the different simulations have generally good agreement in the likely location during gaps, while in others they diverge. Multiple imputation can be particularly powerful if we want to incorporate environmental variables, as spatially explicit variables can be extracted for each simulated track to sample the most likely conditions encountered by the animal.

5 Filtering error

5.1 Gentoo movement data

We have the data from 5 gentoo penguins. We load it and change the ID to a factor and the datetime into a datetime object

gentoo_tracks <- read.csv(
  here("D1-data-prep-ssm", "data", "tracks_argos.csv")) %>%
  mutate(ID = factor(ID),
         datetime = ymd_hms(datetime))
  
gentoo_tracks <- gentoo_tracks %>% subset(ID != "170895" & ID != "170936")
head(gentoo_tracks)

##          ID            datetime Longitude  Latitude Qual
## 1218 170902 2018-05-24 02:04:31 -59.19942 -51.41017    0
## 1219 170902 2018-05-24 02:04:31 -59.19942 -51.41017    0
## 1220 170902 2018-05-24 02:04:31 -59.19942 -51.41017    0
## 1221 170902 2018-05-24 02:04:31 -59.19942 -51.41017    0
## 1222 170902 2018-05-24 10:09:38 -59.23020 -51.41526    2
## 1223 170902 2018-05-24 10:09:38 -59.23020 -51.41526    2

The columns/variables are: - ID: Individual identifier - time: time of location - x: longitude - y: latiture - loc_class: all GPS

We want to remove any row that is missing latitude or longitude data.

gentoo_tracks <- gentoo_tracks %>% 
  filter(!is.na(Longitude) & !is.na(Latitude),
         # remove identical records
         !(datetime == lag(datetime) & 
             Longitude == lag(Longitude) & 
             Latitude == lag(Latitude) & 
             Qual == lag(Qual)))

Here, the data are already sorted in order of ID and datetime, but this is not always the case.

Let’s make a spatial object and plot the data.

gentoo_sf <- gentoo_tracks %>%
  st_as_sf(coords = c("Longitude", "Latitude")) %>% # converts to an sf object
  st_set_crs(4326)

ggplot() +
  geom_sf(data = gentoo_sf, aes(col = ID)) 

Let’s look at the data gaps as above.

# Calculate time difference between locations
gentoo_tracks <- gentoo_tracks %>%
  mutate(dt = ifelse(ID == lead(ID), # If next data row is same individual
                     # calculate time difference
                     difftime(lead(datetime), datetime, units = "mins"), 
                     NA))

# Visualise time differences (all and zoomed)
par(mfrow = c(1, 2))
hist(gentoo_tracks$dt, 1000, main = NA, xlab = "Time difference (min)")
hist(gentoo_tracks$dt, 1000, main = NA, xlab = "Time difference (min)",
     xlim = c(0,2*60))

ggplot() +
  geom_line(data = gentoo_tracks, aes(x = datetime, y = dt/60)) +
  facet_wrap(~ID) +
  ylab("Time interval (hour)") +
  xlab("Date")

# identify the most frequent dt
gentoo_tracks %>% 
  {table(round(.$dt))} %>% 
  sort(decreasing = TRUE) %>% 
  head(20)
## 
##   0   1   2   6 239   3   4  12  17   5  96  97  11 103 242 487 722   7  10  13 
##  37  25  22  18  15  14  12  11  11  10  10  10   9   9   8   8   8   7   7   7

We see that the most frequent time gap is \(0\) min, followed by \(2\), \(1\), and \(6\) min. But we also see a large variability.

We can now use the function p_na (in the script utility_functions.R) to look at the proportion of NAs we would get with 2hr (120 min), 4 hr (240 min), and 6 hr (360), 12 hr (720) resolution.

dt_s <- c(2, 4, 6, 12) # hours
# summarise track dt
gentoo_tracks %>% 
  group_by(ID) %>% 
  summarise(p_NA_1 = p_na(datetime, dt_s[1]*60),    
            p_NA_2 = p_na(datetime, dt_s[2]*60),    
            p_NA_3 = p_na(datetime, dt_s[3]*60),  
            p_NA_4 = p_na(datetime, dt_s[4]*60)) %>% 
  # return formatted table
  kable(digits = 3, 
        col.names = c("Gentoo ID", 
                      paste(dt_s, "h"))) %>%
  kable_styling() %>% 
  add_header_above(c("", "Resolution" = 4))

	Resolution
Gentoo ID	2 h	4 h	6 h	12 h
170902	0.705	0.530	0.406	0.082
170905	0.730	0.532	0.421	0.103
170914	0.731	0.546	0.427	0.138

The proportions of NAs estimated here is very approximate. Worth looking further into this as its common for Argos data to have many locations in a row (i.e., clumps).

We may also want to split first. From the time series, it looks like the biggest gaps are in the order of 20 hr or more. So we are going to split at 20 hr. The function split_at_gap require the column name ID and time. So we reformat before.

gentoo_tracks <- gentoo_tracks %>% rename(time = datetime)
# Use function from utility_function.R to split data at gaps
gentoo_split <- split_at_gap(data = gentoo_tracks, 
                           max_gap = 20*60, 
                           shortest_track = 20*60)
head(gentoo_split)

##         ID                time Longitude  Latitude Qual          dt ID_old
## 1 170902-1 2018-05-24 10:09:38 -59.23020 -51.41526    2   1.4833333 170902
## 2 170902-1 2018-05-24 10:11:07 -59.23301 -51.41519    B 198.7000000 170902
## 3 170902-1 2018-05-24 13:29:49 -59.19302 -51.38425    B 762.7166667 170902
## 4 170902-1 2018-05-25 02:12:32 -59.23865 -51.41388    0 482.0333333 170902
## 5 170902-1 2018-05-25 10:14:34 -59.23944 -51.42496    A  76.7833333 170902
## 6 170902-1 2018-05-25 11:31:21 -59.20930 -51.42498    B   0.9833333 170902

unique(gentoo_split$ID)

##  [1] "170902-1"  "170902-2"  "170902-3"  "170902-5"  "170902-6"  "170902-7" 
##  [7] "170902-8"  "170905-1"  "170905-2"  "170905-5"  "170905-6"  "170905-7" 
## [13] "170905-8"  "170905-9"  "170905-10" "170914-1"  "170914-3"  "170914-4" 
## [19] "170914-5"  "170914-6"  "170914-8"  "170914-9"  "170914-10" "170914-11"
## [25] "170914-13" "170914-15" "170914-16" "170914-17"

dt_s <- c(2, 4, 6, 12) # hours
# summarise track dt
gentoo_split %>% 
  group_by(ID) %>% 
  summarise(p_NA_1 = p_na(time, dt_s[1]*60),    
            p_NA_2 = p_na(time, dt_s[2]*60),    
            p_NA_3 = p_na(time, dt_s[3]*60),  
            p_NA_4 = p_na(time, dt_s[4]*60)) %>% 
  # return formatted table
  kable(digits = 3, 
        col.names = c("Gentoo ID", 
                      paste(dt_s, "h"))) %>%
  kable_styling() %>% 
  add_header_above(c("", "Resolution" = 4))

	Resolution
Gentoo ID	2 h	4 h	6 h	12 h
170902-1	0.600	0.500	0.250	0.000
170902-2	0.704	0.529	0.382	0.029
170902-3	0.717	0.500	0.400	0.056
170902-5	0.669	0.485	0.351	0.052
170902-6	0.659	0.467	0.355	0.000
170902-7	0.636	0.435	0.250	0.000
170902-8	0.593	0.357	0.263	0.000
170905-1	0.559	0.389	0.250	0.000
170905-10	0.706	0.494	0.333	0.034
170905-2	0.697	0.505	0.409	0.059
170905-5	0.700	0.455	0.286	0.000
170905-6	0.690	0.467	0.364	0.000
170905-7	0.689	0.381	0.267	0.000
170905-8	0.677	0.447	0.344	0.059
170905-9	0.703	0.432	0.292	0.000
170914-1	0.700	0.429	0.300	0.000
170914-10	0.731	0.432	0.269	0.000
170914-11	0.650	0.400	0.375	0.000
170914-13	0.625	0.267	0.200	0.000
170914-15	0.729	0.556	0.417	0.077
170914-16	0.531	0.320	0.154	0.000
170914-17	0.717	0.458	0.375	0.000
170914-3	0.455	0.167	0.200	0.000
170914-4	0.654	0.462	0.333	0.000
170914-5	0.741	0.500	0.368	0.000
170914-6	0.583	0.429	0.200	0.000
170914-8	0.737	0.536	0.368	0.000
170914-9	0.688	0.375	0.333	0.000

The results change when splitting, but it’s still hard to find a perfect resolution.

5.2 Fit state-space model

Filter out error, but return at same irregular time interval using a correlated random walk

gentoo_mpi <- fit_ssm(gentoo_tracks,
                      id = "ID",
                      coord = c("Longitude", "Latitude"),
                      date = "time",
                      lc = "Qual",
                      model = "mp")

## fitting mp SSM to 3 tracks...
##  pars:   0.78793 1.45317 0 -0.20649 0 0 0       pars:   0.58693 1.10913 -0.06442 -0.21683 -0.68041 -0.61155 -0.00536       pars:   -0.01607 0.07702 -0.25768 -0.24785 -2.72165 -2.44619 -0.02144       pars:   1.15451 -0.62144 -1.76269 -0.47685 -0.02191 -4.40014 0.82588       pars:   0.12016 0.12577 -0.34236 -0.26059 -2.34668 -2.40817 0.0379       pars:   0.10272 0.21628 -0.6613 -0.3158 -2.13094 -2.39738 0.15584       pars:   -0.01744 0.16348 -0.99362 -0.40382 -2.29233 -2.45403 0.2516       pars:   0.10232 0.21039 -1.27296 -0.48171 -2.06454 -2.34889 0.35284       pars:   0.03201 0.06932 -1.49576 -0.6253 -2.30079 -2.45309 0.46268       pars:   0.10052 0.3739 -1.61469 -0.82414 -2.27997 -2.38691 0.59641       pars:   0.14314 0.16671 -1.75348 -0.95266 -2.06767 -2.56331 0.72726       pars:   -0.05183 0.11061 -1.84963 -1.04434 -2.18817 -2.27515 0.85945       pars:   0.1886 0.31714 -2.00288 -1.10365 -2.21214 -2.46379 0.96056       pars:   -0.09481 0.25872 -2.11849 -1.28867 -2.12676 -2.63587 1.03544       pars:   0.08361 0.24179 -2.06631 -1.12723 -2.22119 -2.48542 0.98443       pars:   0.08729 0.30211 -2.1481 -1.13715 -2.13752 -2.41613 0.99689       pars:   0.09021 0.20881 -2.24786 -1.17144 -2.18392 -2.41103 1.01628       pars:   0.0412 0.32168 -2.32098 -1.18689 -2.17151 -2.4409 1.04127       pars:   0.07656 0.27552 -2.42925 -1.21263 -2.12499 -2.46831 1.10223       pars:   0.09825 0.2833 -2.53925 -1.22897 -2.2019 -2.41179 1.12197       pars:   0.04101 0.26608 -2.62224 -1.30868 -2.15866 -2.45203 1.07619       pars:   0.11457 0.34924 -2.68725 -1.28006 -2.10994 -2.42122 1.11747       pars:   0.06397 0.29204 -2.64253 -1.29974 -2.14345 -2.44241 1.08908       pars:   0.06022 0.27825 -2.66999 -1.29675 -2.16405 -2.45986 1.11114       pars:   0.05865 0.29043 -2.71043 -1.25886 -2.13422 -2.44542 1.17755       pars:   0.06895 0.28918 -2.7808 -1.31665 -2.15121 -2.43979 1.17395       pars:   0.06895 0.28918 -2.7808 -1.31665 -2.15121 -2.43979 1.17395      
##  pars:   1.52139 1.9791 0 0.80904 0 0 0       pars:   1.22717 1.5513 -0.08327 0.86376 -0.63699 -0.56094 0.00926       pars:   0.34454 0.2679 -0.33307 1.02789 -2.54797 -2.24376 0.03705       pars:   -0.20967 -0.45248 -0.61837 1.14349 -2.70732 -2.85617 0.04683       pars:   0.95063 1.16564 -0.18626 0.91759 -1.03546 -1.00269 0.01649       pars:   0.36826 0.26103 -0.59304 1.0546 -1.38774 -1.95395 0.02235       pars:   -0.38394 0.50359 -1.36545 1.11435 -0.47491 -1.4643 -0.24769       pars:   -0.51669 0.43549 -1.7905 1.20048 -1.06486 -1.61772 -0.35142       pars:   -0.01875 0.22151 -1.94238 1.21716 -0.70944 -2.00338 -0.3276       pars:   -0.37169 0.35202 -1.8222 1.20217 -0.84305 -1.73104 -0.33508       pars:   -0.22751 0.38614 -1.98148 1.25256 -1.03488 -1.7818 -0.28966       pars:   -0.09649 0.47497 -2.19359 1.31861 -0.92257 -1.73806 -0.24104       pars:   -0.1161 0.47779 -2.36395 1.5281 -0.97105 -1.85708 -0.2796       pars:   -0.16764 0.39241 -2.46766 1.64813 -1.03234 -1.62815 -0.28467       pars:   -0.13836 0.48804 -2.36717 1.53196 -0.98449 -1.78459 -0.28143       pars:   -0.13752 0.41699 -2.38319 1.55353 -0.96956 -1.78476 -0.28996       pars:   -0.15527 0.45028 -2.38277 1.5529 -0.97392 -1.77615 -0.28993       pars:   -0.15527 0.45028 -2.38277 1.5529 -0.97392 -1.77615 -0.28993      
##  pars:   1.5194 2.01369 0 -1.63902 0 0 0       pars:   1.17407 1.55397 -0.07777 -1.73645 -0.59921 -0.54274 -0.01566       pars:   0.13807 0.17483 -0.31109 -2.02872 -2.39685 -2.17097 -0.06262       pars:   1.91168 -0.88118 -0.90157 -1.45951 0.33805 -3.91506 0.67559       pars:   0.71613 0.32201 -0.3612 -1.8531 -1.45622 -2.08779 0.08339       pars:   0.52107 -0.24337 -1.06108 -2.0901 -1.18801 -1.54091 0.24909       pars:   0.63474 0.08612 -0.65321 -1.95198 -1.34431 -1.85962 0.15252       pars:   0.67551 0.20428 -0.50693 -1.90245 -1.40037 -1.97392 0.11789       pars:   0.59217 0.25545 -0.54589 -1.88632 -1.367 -2.01213 0.13196       pars:   0.33848 0.32264 -0.74762 -1.87169 -1.23247 -2.04687 0.19425       pars:   -0.04116 0.32153 -1.14495 -1.88934 -0.99334 -2.00521 0.30978       pars:   -0.4679 0.04848 -1.40239 -1.95427 -1.59389 -2.6754 0.34333       pars:   -0.15947 0.24583 -1.21632 -1.90734 -1.15983 -2.19101 0.31908       pars:   -0.23328 0.21674 -1.42964 -1.89946 -1.06154 -2.35201 0.32381       pars:   -0.21829 0.26265 -1.72036 -1.91211 -1.04279 -2.3565 0.31433       pars:   -0.29254 0.408 -1.87497 -1.85051 -1.20041 -2.44585 0.29702       pars:   -0.28425 0.26038 -1.951 -1.65257 -1.1145 -2.44949 0.1819       pars:   -0.29715 0.32419 -2.00685 -1.66196 -1.1419 -2.16923 0.15494       pars:   -0.2965 0.29226 -1.98077 -1.66482 -1.12188 -2.35729 0.1778       pars:   -0.28893 0.28537 -2.07396 -1.69201 -1.10217 -2.38119 0.19475       pars:   -0.42606 0.45421 -2.26166 -1.82373 -1.04706 -2.42058 0.30284       pars:   -0.11068 0.02842 -2.09264 -1.80935 -1.12488 -2.45073 -0.04007       pars:   -0.29164 0.25749 -2.06799 -1.70557 -1.12744 -2.38858 0.1896       pars:   -0.28954 0.27907 -2.07261 -1.69507 -1.10788 -2.38286 0.19359       pars:   -0.28976 0.28424 -2.07756 -1.69791 -1.1118 -2.38644 0.19434       pars:   -0.29101 0.2795 -2.07475 -1.70485 -1.11329 -2.38681 0.19196       pars:   -0.28985 0.28204 -2.07768 -1.71251 -1.11436 -2.39007 0.19083       pars:   -0.28985 0.28204 -2.07768 -1.71251 -1.11436 -2.39007 0.19083

head(gentoo_mpi)

## # A tibble: 3 × 5
##   id     ssm           converged pdHess pmodel
##   <chr>  <named list>  <lgl>     <lgl>  <chr> 
## 1 170902 <mp_ssm [15]> TRUE      TRUE   mp    
## 2 170905 <mp_ssm [15]> TRUE      TRUE   mp    
## 3 170914 <mp_ssm [15]> TRUE      TRUE   mp

gentoo_mpi

## # A tibble: 3 × 5
##   id     ssm           converged pdHess pmodel
##   <chr>  <named list>  <lgl>     <lgl>  <chr> 
## 1 170902 <mp_ssm [15]> TRUE      TRUE   mp    
## 2 170905 <mp_ssm [15]> TRUE      TRUE   mp    
## 3 170914 <mp_ssm [15]> TRUE      TRUE   mp

summary(gentoo_mpi)

##  Animal id Model Time n.obs n.filt n.fit n.pred n.rr converged   AICc
##     170902    mp    .   404     31   373      0    .      TRUE 4879.6
##     170905    mp    .   305     26   279      0    .      TRUE 4325.6
##     170914    mp    .   350     55   295      0    .      TRUE 4223.7
## 
## --------------
## 170902 
## --------------
##  Parameter Estimate Std.Err
##      rho_p  -0.5772  0.0753
##    sigma_x   1.0714  0.0827
##    sigma_y   1.3353   0.086
##      rho_o   0.5277  0.0795
##      tau_x   0.1163  0.0076
##      tau_y   0.0872  0.0066
##    sigma_g    0.062  0.0145
## 
## --------------
## 170905 
## --------------
##  Parameter Estimate Std.Err
##      rho_p   0.6507  0.0949
##    sigma_x   0.8562  0.1154
##    sigma_y   1.5687  0.1337
##      rho_o   -0.144  0.1157
##      tau_x   0.3776  0.0257
##      tau_y   0.1693  0.0126
##    sigma_g   0.0923  0.0248
## 
## --------------
## 170914 
## --------------
##  Parameter Estimate Std.Err
##      rho_p  -0.6943  0.0789
##    sigma_x   0.7484  0.0769
##    sigma_y   1.3258  0.1017
##      rho_o   0.0951  0.1299
##      tau_x   0.3281    0.02
##      tau_y   0.0916  0.0099
##    sigma_g   0.1252  0.0292

Done individual. “Currently, this model can only be fit to individuals separately.”

plot(gentoo_mpi, what ="fitted", type=2, ask = FALSE)

## $`170902`

## 
## $`170905`

## 
## $`170914`

plot(gentoo_mpi, what ="fitted", type=1, ask = FALSE)

## $`170902`

## 
## $`170905`

## 
## $`170914`

Observations that fail the prefilter stage are displayed with x. Rug are observation.

Map the behavioural persistance model.

aniMotum::map(gentoo_mpi, what = "fitted")

## using map scale: 10

“To aid visualisation of high and low move persistence regions along each track, the γt estimates are normalised separately for each track to span the interval 0,1. This can be turned off with normalise = FALSE, or applied to the group of tracks with group = TRUE.” “When fitting to a collection of individuals, the normalisation can be applied separately to each individual or collectively via normalise = TRUE, group = TRUE. In the latter case, the relative magnitudes of move persistence are preserved across individuals.”

aniMotum::map(gentoo_mpi, what = "fitted", normalise = TRUE, group = TRUE)

## using map scale: 10

All individual mapped.

TO be able to analyse the data further or to be able to plot it yourself, you can graph the data and make it into a

gentoo_mpi_t <- grab(gentoo_mpi)
head(gentoo_mpi_t)

## # A tibble: 6 × 11
##   id     date                  lon   lat      x      y  x.se  y.se logit_g
##   <chr>  <dttm>              <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1 170902 2018-05-24 10:09:38 -59.2 -51.4 -6593. -6662. 0.179 0.112   0.652
## 2 170902 2018-05-24 10:11:07 -59.2 -51.4 -6593. -6662. 0.183 0.123   0.652
## 3 170902 2018-05-25 02:12:32 -59.2 -51.4 -6594. -6661. 1.49  1.26    0.652
## 4 170902 2018-05-25 10:14:34 -59.2 -51.4 -6593. -6663. 1.64  1.31    0.652
## 5 170902 2018-05-25 11:31:21 -59.2 -51.4 -6592. -6663. 1.57  1.23    0.652
## 6 170902 2018-05-25 11:32:20 -59.2 -51.4 -6592. -6663. 1.56  1.22    0.652
## # ℹ 2 more variables: logit_g.se <dbl>, g <dbl>

gentoo_mpi_sf <- gentoo_mpi_t %>%
  st_as_sf(coords = c("lon", "lat")) %>% # converts to an sf object
  st_set_crs(4326)

ggplot() +
  geom_spatial_path(data =  gentoo_mpi_t, 
                    aes(x = lon, y = lat, lty = id), 
                    crs = 4326) +
  geom_sf(data = gentoo_mpi_sf, aes(col=g)) 

Filter and regularize

Given how

gentoo_mpr <- fit_ssm(gentoo_tracks,
                      id = "ID",
                      coord = c("Longitude", "Latitude"),
                      date = "time",
                      lc = "Qual",
                      model = "mp",
                      time.step = 10)

## fitting mp SSM to 3 tracks...
##  pars:   0.62266 1.18105 0 -0.19648 0 0 0       pars:   0.40072 0.85812 -0.05583 -0.20596 -0.68409 -0.61258 -0.00433       pars:   -0.26509 -0.11069 -0.22333 -0.23439 -2.73636 -2.45031 -0.0173       pars:   1.66152 -0.53761 -1.64874 -0.53655 -0.31922 -4.39194 0.59123       pars:   -0.00934 -0.02063 -0.32579 -0.25704 -2.25554 -2.37832 0.03754       pars:   -0.06208 0.22025 -0.78645 -0.38201 -2.16816 -2.33435 0.1979       pars:   -0.0628 -0.07047 -1.11526 -0.49012 -2.25327 -2.64889 0.31368       pars:   -0.06244 0.07528 -0.95041 -0.43592 -2.2106 -2.4912 0.25564       pars:   -0.03628 0.21809 -1.14952 -0.50854 -2.19744 -2.40362 0.33749       pars:   0.04791 0.05682 -1.58154 -0.67343 -2.19209 -2.59236 0.53647       pars:   0.07379 0.23342 -1.83444 -0.98551 -2.20683 -2.39909 0.83878       pars:   -0.12843 0.31067 -1.8754 -1.29412 -2.09078 -2.63992 1.16578       pars:   0.00157 0.1959 -1.88118 -1.07854 -2.20626 -2.53175 0.93998       pars:   0.07827 0.17202 -2.04355 -1.12799 -2.2116 -2.43748 0.97359       pars:   -0.00972 0.27791 -2.16376 -1.16648 -2.26095 -2.50996 1.02585       pars:   0.04996 0.22878 -2.19071 -1.1684 -2.19311 -2.50773 1.02966       pars:   0.03839 0.23117 -2.26602 -1.18935 -2.23786 -2.45632 1.05045       pars:   0.00965 0.27997 -2.44934 -1.20373 -2.15245 -2.48265 1.06772       pars:   0.06157 0.15147 -2.5667 -1.22517 -2.22925 -2.54985 1.10579       pars:   0.02337 0.21433 -2.47233 -1.21437 -2.1954 -2.50218 1.08024       pars:   0.07355 0.26111 -2.51807 -1.22296 -2.21188 -2.48567 1.09055       pars:   0.04513 0.22894 -2.58406 -1.23655 -2.19009 -2.47956 1.11565       pars:   0.04136 0.24748 -2.64588 -1.26153 -2.23568 -2.46387 1.13551       pars:   0.04326 0.23813 -2.61469 -1.24893 -2.21268 -2.47179 1.12549       pars:   0.06107 0.23239 -2.64561 -1.25943 -2.19777 -2.48515 1.13089       pars:   0.04826 0.24593 -2.6792 -1.27019 -2.20099 -2.47134 1.13825       pars:   0.04826 0.24593 -2.6792 -1.27019 -2.20099 -2.47134 1.13825      
##  pars:   1.36081 1.63992 0 0.6902 0 0 0       pars:   1.06008 1.27447 -0.06057 0.73308 -0.62385 -0.61726 0.01731       pars:   0.1579 0.17812 -0.24228 0.86171 -2.4954 -2.46905 0.06924       pars:   -0.27876 -0.47656 -0.62302 0.95676 -2.66894 -3.18472 0.04546       pars:   0.81864 0.9587 -0.162 0.77342 -0.99265 -1.08026 0.02239       pars:   0.32328 0.24274 -0.40601 0.85797 -1.48196 -2.08859 0.00094       pars:   -0.28254 -0.00286 -0.98215 0.90527 -0.41391 -1.74101 -0.24039       pars:   -1.01896 0.3265 -1.71209 1.03909 -1.16531 -1.20289 -0.40038       pars:   -0.4203 0.13449 -1.07615 0.9396 -1.02379 -1.70181 -0.23038       pars:   -0.3615 0.32269 -1.60685 1.08032 -0.85238 -1.912 -0.14346       pars:   -0.36668 0.13202 -2.12218 1.25947 -0.96847 -1.64258 -0.10689       pars:   -0.48605 0.17035 -2.02328 1.29981 -1.00083 -1.91744 -0.14737       pars:   -0.61883 0.36448 -1.9897 1.34938 -0.85541 -1.77169 -0.21113       pars:   -0.46521 0.22627 -2.01846 1.2992 -0.93141 -1.80659 -0.15131       pars:   -0.37757 0.23535 -1.94051 1.36651 -0.97537 -1.81767 -0.16812       pars:   -0.43405 0.25796 -1.89013 1.45719 -0.92685 -1.81217 -0.23025       pars:   -0.42169 0.16415 -1.88251 1.54861 -0.91838 -1.77763 -0.27579       pars:   -0.43405 0.25796 -1.89013 1.45719 -0.92685 -1.81217 -0.23025      
##  pars:   1.27091 1.67052 0 -1.25101 0 0 0       pars:   0.95858 1.26026 -0.04637 -1.32967 -0.61765 -0.58669 -0.01041       pars:   0.02161 0.0295 -0.18548 -1.56567 -2.47059 -2.34677 -0.04165       pars:   -0.38873 -0.57955 -0.4983 -1.68877 -2.71622 -3.04841 -0.06097       pars:   0.59268 0.7606 -0.16911 -1.4272 -1.18758 -1.25525 -0.02414       pars:   -0.20652 -0.43794 -0.69748 -1.63724 -1.62003 -2.71503 0.03827       pars:   -1.37275 0.9647 -1.40339 -1.4443 -0.74364 -2.52437 0.2684       pars:   0.00057 -0.05094 -1.94305 -0.71214 -0.70851 -1.54344 0.37853       pars:   -0.3817 -0.1053 -1.88413 -0.79081 -1.48182 -2.19357 0.40426       pars:   -0.12861 0.21865 -1.59295 -0.98718 -0.82656 -2.82902 0.17062       pars:   -0.31259 0.06343 -1.86652 -0.82821 -1.04917 -2.21545 0.35975       pars:   -0.30032 0.15455 -1.70384 -1.02811 -1.2152 -2.5401 0.2303       pars:   -0.39359 0.28501 -1.7216 -1.31278 -1.22027 -2.25487 0.04018       pars:   -0.32813 0.24952 -1.72286 -1.06629 -1.1592 -2.36795 0.19696       pars:   -0.41505 0.13422 -1.71597 -1.20189 -1.17943 -2.39815 0.12703       pars:   -0.34586 0.18199 -1.79588 -1.37233 -1.1331 -2.41595 0.10415       pars:   -0.33882 0.19862 -1.75447 -1.56842 -1.20232 -2.40147 0.10517       pars:   -0.35911 0.17965 -1.77111 -1.45446 -1.18952 -2.41507 0.10473       pars:   -0.37302 0.16533 -1.71123 -1.50779 -1.13746 -2.39967 0.13555       pars:   -0.33875 0.17761 -1.75381 -1.5841 -1.14406 -2.40433 0.17674       pars:   -0.33875 0.17761 -1.75381 -1.5841 -1.14406 -2.40433 0.17674

head(gentoo_mpr)

## # A tibble: 3 × 5
##   id     ssm           converged pdHess pmodel
##   <chr>  <named list>  <lgl>     <lgl>  <chr> 
## 1 170902 <mp_ssm [15]> TRUE      TRUE   mp    
## 2 170905 <mp_ssm [15]> TRUE      TRUE   mp    
## 3 170914 <mp_ssm [15]> TRUE      TRUE   mp

plot(gentoo_mpr, what ="predicted", type=1, ask=FALSE)

## $`170902`

## 
## $`170905`

## 
## $`170914`

plot(gentoo_mpr, what ="predicted", type=2, ask=FALSE)

## $`170902`

## 
## $`170905`

## 
## $`170914`

aniMotum::map(gentoo_mpr, what = "predicted")

## using map scale: 10

To be able to do multiple imputation, we could simulate from the posterior. Here we do 4, just to be able to vizualise it.

“Posterior simulations from an SSM fit provide a useful characterization of track uncertainty. The sim_post function provides an efficient method for simulating tracks from the joint uncertainty of the estimated locations.”

psim <- sim_post(gentoo_mpr, "p", reps = 4)
gentoo_post <- psim %>% 
  unnest(cols = psims) %>% 
  mutate(rep = factor(rep))

Plots

gentoo_post_sf <- gentoo_post %>%
  st_as_sf(coords = c("lon", "lat")) %>% # converts to an sf object
  st_set_crs(4326)

ggplot() +
  geom_spatial_path(data = gentoo_post, 
                    aes(x = lon, y = lat, colour = rep),
                    lwd = 0.3,
                    crs = 4326) +
  geom_sf(data = gentoo_post_sf, aes(pch = id, colour = rep), size=0.5)

We can see the slightly different paths.

One-step ahead residuals.

res.crw <- osar(gentoo_mpi)
plot(res.crw, type = "ts")

## Warning: Removed 12 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 25 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 31 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 31 rows containing missing values or values outside the scale range
## (`geom_point()`).

plot(res.crw, type = "qq")

plot(res.crw, type = "acf")